Many-body diagrammatic expansion in a Kohn-Sham basis: 
implications for Time-Dependent Density Functional Theory of excited states 
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We formulate diagrammatic rules for many-body perturbation theory which uses Kohn-Sham 
(KS) Green's functions as basic propagators. The diagram technique allows to study the properties 
of the dynamic nonlocal exchange-correlation (xc) kernel f^c- We show that the spatial non- locality 
of fxc is strongly frequency- dependent. In particular, in extended systems the non-locality range 
diverges at the excitation energies. This divergency is related to the discontinuity of the xc potential. 
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Time-dependent density functional theory (TDDFT) 
offers a possibility to extend the powerful density- 
functional formalism to excited states of many-body 
systems [||^. A substantial improvement of excitation 
energies with respect to KS eigenvalues was obtained for 
atoms and molecules |^-^ using a variety of approxima- 
tions for a dynamic xc kernel fxc = 5vxc{Y,t) / 5n{r' ,t') 
{vxc is a xc potential). However in solids the wrong KS 
band gap remains unchanged regardless the approxima- 
tion used, albeit the dielectric function is on average im- 
proved 0. 

This situation, as we show below, reflects an extremely 
nonlocal behavior of fxc at excitation frequencies. The 
non-locality range is as large as the system size and hence 
diverges in extended systems. None of the up-to-date ap- 
proximations account for this behavior as they all employ 
the adiabatic (frequency-independent) xc kernels. 

In this paper we develop a perturbative technique with 
KS Green's functions as the bare propagators. In essence, 
it is a diagrammatic expansion of Sham-Schliiter equa- 
tion which maintains a correct electron density in 
every order of the perturbation theory. We find that at 
resonant frequencies the kernel fxc is proportional to the 
discontinuity of Vxc- This explains the anomalous non- 
locality of fxc, since a constant shift of a potential due 
to an extra particle is felt by another particle anywhere 
in the system. 

In the framework of TDDFT the excitation energies 
are commonly calculated [||-^ from the poles of the lin- 
ear response function x(r, r',w). The latter is related to 
the KS susceptibihty xs(r,r',w) by 



X{^) ^ Xs{^) + Xs{^) Vc + fxc{^^) x(t^), 



(1) 



where Vc = e^/|r — r'| is a bare Coulomb repulsion and 
the kernel fxc enters as an additional dynamic interac- 
tion. 

Alternatively, the poles of x(r, r',^) can be found as 
the eigenvalues of a linearized equation for density matrix 



-ps(ri,r2) j dv\v^{vi,v)-V^{v,V2)\5p{v,v)^{) (2) 

where — Vc + fxd^^), Hs{r) is the KS Hamiltonian, 
and pg{r,r') — rijV'j' (r)^/'j(r') is the equilibrium KS 
density matrix with the KS orbitals ip*{r). Equation 
(§) clearly shows that a correction to the KS excitation 
energies originates from the Hartree-type energy of the 
excitation-induced density fluctuation 6n{r) = 6p{r,r). 
In the KS basis equation (0) takes the form 



kl 



\'^ki)Spki = 0, (3) 



Ef — Ej is a KS excitation energy, 7^ = 
is the difference of the occupation numbers and 
(r) = '0*(r)'0j (r). The ordinary perturbation theory 



where w£- 

n, - 



gives the energy shift AuJij 
as 



in the first order 



,(1) - 



;<i>.,(r)|yc(r,r') + /.c(r,r',c.|.)l'i>..(r')>, (4) 



which is identical to the result of the first-order Laurent 
expansion of x{r, r' ,uj) The perturbation theory eas- 
ily allows to obtain the higher-order corrections as well. 
It is worth noting that equations (H) and (^ and the 
perturbative result (^) are equally valid for finite and for 
extended systems. 

Let us consider the dependence of the first-order cor- 
rection (^) on the size of a system L ~ l/^/^ at fixed av- 
erage density N/V. As ^^(r) contains a normalization 
factor the first term in Eq. (^ is proportional to 

e^/L. This is the Coulomb energy of the density variation 
due to an electron-hole excitation, which is infinitesimally 
small in extended systems. The second term crucially de- 
pends on the non-locality of fxc i-e. on the extension of 
a xc-hole. At w = the non-locality range is about the 
interparticle distance I ~ {V/Ny^^^. This feature is re- 
produced by the popular OEP approximation j^] whereas 
in adiabatic LDA |^ fxc is a point interaction. Assum- 
ing that at resonance frequency fxci^ij) has a similar 
non-locality we find that the second term in Eq. (Eh is 
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proportional to \/N and vanishes as (see, however, 
pof). Thus any xc-kernel which is finite and decays at 
infinity does not contribute to LOij in extended systems. 
Using the many-body perturbative approach formulated 
below, we show that a non- vanishing xc correction arises 
due to the divergence of the non- locality scale of fxci^ij)- 




/■ > g- > h. 
FIG. 2. First (a) and second (b-i) -order diagrams for xc 
potential. 



FIG. 1. First-order corrections to the Green's function. 

We employ the Matsubara formalism at nonzero temper- 
ature T which enables us to obtain any physical retarded 
function through analytic continuation [Tl] ]. Assuming 
that the xc-potential Vxc{^) and the ground state density 
are known, we represent the Hamiltonian of a system as 
a sum of the KS Hamiltonian Hs and the perturbation 
V ^ where 



fc=i 



(5) 



In Eq. (H) W' is a two-particle interaction with the 
Hartree part being subtracted, n(r) is a density operator 
and we assume that the xc-potential v^c can be expanded 
in a power series Vxc — X^fcLi ^^c\ where vi'^c ~ e^*"'. 

Following the standard procedure we define the 
Green's function 



G{X,X') = -{Tr^s{X)^tiX')a)/{a), 



(6) 



where X = (r, r) (r is an imaginary time), ^5(X) is a 
field operator in a KS interaction representation and a 
is a Matsubara S-matrix |ll[] which corresponds to the 
perturbation V of Eq. ^ . The angular brackets in Eq. 
denote averaging over the KS equilibrium state. 
A formal graphical expansion of G{X,X') contains 
along with the pair interaction terms the diagrams re- 

(k) 

lated to the scattering by the "external" potentials Vxc ■ 
To achieve a closed scheme one needs a complemen- 
tary graphical representation of Vxci'fc), which can be 
obtained from the condition of density conservation . 
The Green's function (H) can be written in the form 
G = Gs + J2T=i G'^^^ ' where Gs is a KS Green's function 
and G^'^^ is a fc-th order correction. As the KS system 
possesses an exact density, the variation of the density 
due to the interaction must vanish. Applying this 
requirement to every order we arrive at the system of 
equations: 



(7) 



where ujn = 7rT(2n + 1). This system is equivalent to 
well-known Sham-Schliiter equation |^ (see also A 
successive solution of Eq. (0) allows to construct wic'' 
for every k. For example, the first-order correction G^-^"^ 
is presented in Fig. la, where the solid line is the KS 
Green's function and dashed line is the Coulomb inter- 
action. Substituting G^^^ into Eq. (Q) we get vi]} as 
shown in Fig. 2a, where the wiggled line stands for the 
inverse KS response function X5^(r,r'). The final first- 
order correction to the Green's function is shown in Fig. 
lb. Note that vi]} (Fig. 2a) exactly corresponds to the 
x-only OEP Vx- Given vi^J we solve Eq. (R) for A: = 2 

(2) 

and obtain eight graphs for VxJ (Fig. 2b-i). From the 
further iterations we deduce the following diagrammatic 
rules for Vxc in arbitrary order: 

i. Draw all graphs for density according to the usual 
rules [0 and attach wiggled lines to external point of 
each graph. 

ii. Whenever it is possible separate the graphs into two 
parts by cutting two fermionic lines. Join the external 
fermionic lines of these parts and connect them by the 
wiggled line. Do not cut lines attached to the external 
wiggled line. 

in. If a new graph is separable, repeat ii. 

iv. If there are several possible cross sections repeat ii 

and Hi for every cross-section. 

V. Leave only nonequivalent graphs. 

For example, graphs a-e in Fig. 2 appear according to i. 

Graph f is obtained from d applying ii, whereas graphs g- 

i originate from the graph e by the successive application 

of ii-v. 

Given a diagrammatic representation for Vxc we can 
easily construct a graphical expansion for any quantity 
e.g. for one particle Green's function, a response func- 
tion or an energy. We find that the series for the energy 
coincides with the energy expansion obtained in a dif- 
ferent context in Ref. [|l3| (see also Ref. |lj]). As the 
diagrammatic expansion is derived maintaining the ex- 
act density in every order, the series for Vxc is in fact 
a graphical representation of the Gorling-Levy perturba- 
tion theory (GLPT) An obvious advantage of the 
graphical method is a possibility to construct vi'c' for ev- 
ery fc in a transparent form. 
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An important feature of the KS-based diagram tech- 
nique is that every irreducible self-energy insertion 
S(r,r') is accompanied by the local counter term with 
the opposite sign which has the structure of the aver- 
age of E: {Gs'SGs){GsGs)~'^ ■ This term guarantees 
the density conservation and locally reduces the effective 
field S. The first-order correction (Fig. lb) gives an ex- 
ample of this compensation. It is interesting to note that 
also the standard diagram technique can be reformulated 
in a similar fashion. One has to explicitly introduce the 
correction to the chemical potential to compensate the 
change of the total number of particles (i.e. the averaged 
density) in every order of the perturbation theory. This 
leads to similar, but spatially-averaged counter terms. 
However the local compensation in the KS-based tech- 
nique is obviously more efficient. This means that KS 
particles are much closer to the true quasiparticles than 
bare electrons. 

Our graphical method has an obvious connection to the 
GW approximation [|l6|. Let us collect in every order of 
the perturbation theory only the bubble diagrams (e.g. 
the graph in Fig. 2c) and sum them up. The correspond- 
ing correction to the Green's function is still given by 
Fig. lb, but with the RPA-screened interaction. While 
the first graph in Fig. lb is exactly the GW self-energy, 
the second one deviates from the common GW prescrip- 
tion. Instead of subtraction of the whole v^c, one has to 
use an approximate Vxci (even if the exact Vxc is known!). 
This is necessary for the purpose of internal consistency 
and facilitates the density conservation. From this point 
of view it is clear that KS eigenvalues should describe 
well quasiparticles in metals (e.g. the shape of the Fermi 
surface), but not in semiconductors. Indeed, for a short- 
ranged screened interaction in metals the first (nonlocal) 
and the second (local) term in Fig. 2b almost cancel 
(they would cancel exactly for a point interaction). Con- 
versely, in insulators there is no pronounced cancellation 
since the interaction is long-ranged. As a result, the cor- 
rection to the KS energies gets larger the larger the gap 
is. 

To study the electron-hole excitations one has to con- 
sider the linear response function x(r, i"', w„). An integral 
equation for this function 



X(W«) = X(Wti) + x(^n)^CX(w«), 



(8) 



contains irreducible polarization operator x(a;„). We 
split x[i^n) into two parts: x('^«) = Xs{i^n) + n(w„), 
where n(ti;„) includes all (self-energy and vertex) correc- 
tions to the irreducible susceptibility x(w„). The first- 
order corrections to xi^n) are shown in Fig. 3, where 
the first four terms correspond to the first-order correc- 
tion n^^)(cj„). Thus the total response function in the 
first order is 

X{uJn) ^XsM+^^^^\^n)+Xs{i^n)VcXs{uJn)- (9) 
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FIG. 3. First-order corrections to the response function. 

The graphs in Fig. 3 display the physical meaning of all 
corrections to the electron-hole excitation energy. The 
first two diagrams as well as the third and the fourth 
graphs are the self-energy corrections for an electron and 
a hole respectively. The fifth graph is the electron-hole 
Coulomb interaction, while the last term is the Coulomb 
energy of the excitation-induced charge density. Perform- 
ing summation over frequencies and an analytic continu- 
ation in the first five diagrams in Fig. 3 we obtain II'^"'^^ (uj) 
close to the excitation frequency ujf, : 



nW(^)=^,^.$,^.(r)$,*,.(r')-^^ 



(10) 



A, 



{'iI;,{v)\vf{v, r') - Vx{r)5{v - r')|V.(r')) 
(^,(r)|«;^(r,r') - Vx{v)5{v - r')|V,(r')) 
(<i>.,(r)|Vb(r,r')|<i>,,(r')). (11) 



In Eq. ( |ll| ) vp{y,y') is the Fock nonlocal potential and 
Vx (r) is the local exchange potential or OEP. Apparently 
7ijAy is a shift of the excitation energy due to the first 
five graphs in Fig. 3. The last graph gives a correction 
which is the first term in Eq. (Q). The total correc- 
tion Awlj' — jij{Aij + {^ij\Vc\^ij)) coincides with the 
first-order GLPT With the increase of a sys- 

tem volume the last term in Eq. (^l]) (the energy of the 
electron-hole interaction) vanishes as e^/L. Conversely, 
the first two terms remain finite (of the order of e'^/l) 
and describe e.g. an exchange shift of the band gap of a 
semiconductor. 

The above results can be restated in terms of TDDFT. 
In the linear approximation, i.e. after the first iteration 
Eq. (0) yields 

X(Wn) « XS(W„) + Xs{(^n) [Vc + fxc{uJn)] Xs(Wn)- (12) 

A comparison to Eq. (^l]) leads to the following relation 
fxc{co)^Xs{io)-'U^'Hio)xs{io)-'. (13) 

As the linear approximation vi]} is simply the x-only 
OEP (Fig. 2a), a linear expression Eq. (jl^) should pro- 
vide a dynamical OEP kernel f^^^{uj). Indeed, a direct 
functional differentiation of the graph Fig. 2a versus den- 
sity yields Eq. (|l|). 

At any non-resonant frequency, e.g. in statics lu = 0, 
the non-locality range of fxc Eq. (^) is about an in- 
terparticle distance I. The correction to the excitation 
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energies Eq. (W) depends, however, on the kernel at 
resonance fxd^j)- Let us consider a spatial extension 
of this kernel using Eq. (|l3|). To calculate xs(w)^^ 
at resonance cufj one has to keep both a singular and 
a regular parts of xsi^)- Following Ref. |Q we write 
XS = J^j<^>^jir)<^>^jir')/{uJ - u;fj) + Xr{r,r'), where Xr is a 
regular part. Substituting n(i)(w) (|lO|) to Eq. @ and 
performing calculations we arrive at the following result 

. /rfridr2x7Hi-.i-i)^»j(ri)^rj(r2)x,T'(r2,rO 



[/ dridr2$^(r2)xr'(r2,ri)$y (ri) 



(14) 



It is seen from Eq. (|lj) that a spacial scale of 
fxc{Y, r',ujij) is dictated by the functions $ij (r) which ex- 
tend over the whole volume. Thus the non-locality range 
of the resonant f^^ is simply a system size, which makes 
possible a finite xc contribution {^ij\fxc{^ij)\^ij) = ^ij 
in Eq. (^). This result was considered in Ref. as an 
indication of an equivalence of GLPT and TDDFT. We 
emphasize that this equivalence holds only for a dynamic 
xc-kernel at resonant frequency and is not fulfilled by any 
static approximation. 

It is straightforward to construct a kernel f^^(uj) 
which reproduces the first-order GW result. One has 



to replace n'^^^(a;) in Eq. (13) by the function Uicri^) 
which is defined by the first four graphs in Fig. 3, but 
with the RPA-screened interaction. TDDFT formalism 
with this /^^{lu) exactly reproduces the GW approach. 
For a semiconductor with the band gap Eg the xc-kernel 
at u> = Eg, which is responsible for the band gap correc- 



tion, is given by Eq. (h4 ) with A 



-1— Ajv_i. Here 



A^r-i-i (A^v-i) are the discontinuities of the xc-potential 
upon addition (removal) of a particle. The fundamental 
relation of f^c at resonant frequency to the discontinuity 
of Vxc has a simple physical interpretation. The jump 
of Vxc signifies a constant shift of a potential throughout 
the system due to addition of one particle. It means that 
another probe particle interacts with the first one any- 
where in a system. This should be interpreted as a xc- 
interaction with a length scale of the size of a system and 
with an amplitude equal to the Vxc discontinuity. Impor- 
tantly, this results holds only at resonant frequency when 
a real creation of an electron-hole pair takes place. The 
arguments above show that not only any static approx- 
imation, but also any LDA-based dynamic approxima- 
tion (including any gradient corrections) for fxc cannot 
provide a consistent results for excitation energies and a 
construction of explicit orbital- and frequency-dependent 
functionals similar to Eq. (|l^) is required. A possible 
alternative is a direct calculation of the irreducible po- 
larization operator using the diagram method outlined 



above. This formally allows to express excitation ener- 
gies as functionals of the KS orbitals and consequently of 
the ground state density. 
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